rm(list = ls())

#setwd("~/Dropbox/Energiewende_Project/All_Replication_Files")
setwd("C:/Users/chris/Dropbox/Energiewende_Project/All_Replication_Files")

#Loading in Packages ----------------------------------------------------------------------------------

library(dplyr)
library(sandwich)
library(lmtest)
library(margins)
library(ggplot2)
library(scales)
library(gridExtra)
library(readr)
library(stargazer)
library(reshape2) 
library(kableExtra)
library(knitr)
library(modelsummary)
library(grid)
options(knitr.table.format = "latex")

####################################################################################################################################
############################################Investigate State Management Covariates####################################################################
####################################################################################################################################
###########Prep Data:####
data <- read_csv("data.coal.cleaned.final.csv", col_types = cols(.default = "?"))

state.man.test.d <- lm(supportintervention ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker, data=data)

#state intervention general: female, younger, poorer, east germany, blue collar workers

state.man.test.econ.d <- lm(supportinterventionecon ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker, data=data)

#state intervention for economy (price and growth): poorer, female, younger, east german

state.man.test.frail.d <- lm(supportinterventionfrail ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker, data=data)

#state intervention for frail (elderly and sick): poorer, female, age doesn't matter, east german, blue collar workers

state.man.test.workers.d <- lm(supportinterventionworkers ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker, data=data)

#state intervention for workers/working-age unemployed: poorer, female, younger, east german doesn't matter, blue collar workers


#table of demographics

table1 <- modelsummary(list("Belief State Intervention" = state.man.test.d, 
                            "Just Economy"=state.man.test.econ.d,
                            "Just Sick and Elderly" =state.man.test.frail.d,
                            "Just Workers" = state.man.test.workers.d),
                       gof_omit = "AIC|BIC|Log.Lik|RMSE",
            coef_rename=c("incomequint2" = "Income Quintile 2",
                        "incomequint3" = "Income Quintile 3",
                        "incomequint4" = "Income Quintile 4",
                        "incomequint5" = "Income Quintile 5",
                        "male" = "Male",
                        "under30" = "Age under 30",
                        "btw30and65" = "Age between 30 and 65",
                        "eastgermany" = "East German",
                        "whitecollarworker" = "White Collar Worker"),
             title="Multivariate Regressions, DV: Belief in State Intervention",
             output="latex"
)

regtable1 <- table1 %>%
  kable_styling(font_size = 7)


save_kable(regtable1, "Tables/state_management_demographic_Covariates.tex")



data$party_2021_2 <- relevel(as.factor(data$party_2021_2), "SPD (Sozialdemokratische Partei Deutschlands)")
#interactions with other attitude variables


#0 means left, 10 means right for left-right variable

state.man.test.d.p <- lm(supportintervention ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker + party_2021_2+ left_right_4 + climate_policyn + climatehumans, data=data)

#state intervention general: female, younger, poorer, east germany, blue collar workers

state.man.test.econ.d.p <- lm(supportinterventionecon ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker+ party_2021_2+ left_right_4 + climate_policyn + climatehumans, data=data)

#state intervention for economy (price and growth): poorer, female, younger, east german

state.man.test.frail.d.p <- lm(supportinterventionfrail ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker+ party_2021_2+ left_right_4 + climate_policyn + climatehumans, data=data)

#state intervention for frail (elderly and sick): poorer, female, age doesn't matter, east german, blue collar workers

state.man.test.workers.d.p <- lm(supportinterventionworkers ~ incomequint2 +incomequint3 + incomequint4 + incomequint5 +  male + under30 + btw30and65 + eastgermany+ whitecollarworker+ party_2021_2 + left_right_4 + climate_policyn + climatehumans, data=data)

table2 <- modelsummary(list("Belief State Intervention" = state.man.test.d.p, 
                  "Just Economy"=state.man.test.econ.d.p,
                  "Just Sick and Elderly" =state.man.test.frail.d.p,
                  "Just Workers" = state.man.test.workers.d.p),
          coef_rename =c("incomequint2" = "Income Quintile 2",
                         "incomequint3" = "Income Quintile 3",
                         "incomequint4" = "Income Quintile 4",
                         "incomequint5" = "Income Quintile 5",
                         "male" = "Male",
                         "under30" = "Age under 30",
                         "btw30and65" = "Age between 30 and 65",
                         "eastgermany" = "East German",
                         "whitecollarworker" = "White Collar Worker",
                         "party_2021_2AfD (Alternative für Deutschland)" = "AfD Voter",
                         "party_2021_2Andere Partei" = "Other Party",
                         "party_2021_2CDU/CSU (Christlich Demokratische Union/Christlich-Soziale Union)" = "CDU/CSU Voter",
                         "party_2021_2FDP (Freie Demokratische Partei)" = "FDP Voter",
                         "party_2021_2GRÜNE (Bündnis 90/Die Grünen)" = "Gruene Voter",
                         "party_2021_2Ich bin nicht wahlberechtigt" = "Not Eligible to Vote",
                         "party_2021_2Ich werde nicht wählen"= "Non-Voter",
                         "party_2021_2LINKE (Die Linke)" = "Left Party Voter",
                         "left_right_4" = "Right-wing (continuous)",
                         "climate_policyn" = "Support Climate Policy",
                         "climatehumans" = "Belief In Man-Made Climate Change"),
           title="Multivariate Regressions, Demographics and Other Beliefs",
          output="latex"
          
)

regtable2 <- table2 %>%
  kable_styling(font_size = 7)

save_kable(regtable2 , "Tables/state_management_demographic_policy_covariates.tex")

                         


#binary regressions with plots 

runreg_state_man_cov <- function(yvar, xvar){
  
  fmla.reg = as.formula(paste0(yvar, '~', xvar))
  
  regout <- lm(fmla.reg,data = data) 
  
  reg.name.coef.se <- c(xvar,regout$coefficients[2],sqrt(diag(vcov(regout))[2]) )
  
  return(reg.name.coef.se)
}

runreg_state_man_table<- function(yvar, xvar){
  
  fmla.reg = as.formula(paste0(yvar, '~', xvar))
  
  regout <- lm(fmla.reg,data = data) 

  return(regout)
}

data$incomequintile <- dplyr::recode(data$incomequintiles,
                              "quint1" = 1,
                              "quint2" = 2,
                              "quint3" = 3,
                              "quint4" = 4,
                              "quint5" = 5)

data$fdp <- ifelse(data$party_2021_2=="FDP (Freie Demokratische Partei)",1,0)
data$linke <- ifelse(data$party_2021_2=="LINKE (Die Linke)",1,0)
data$spd <- ifelse(data$party_2021_2=="SPD (Sozialdemokratische Partei Deutschlands)",1,0)
data$cdu <- ifelse(data$party_2021_2=="CDU/CSU (Christlich Demokratische Union/Christlich-Soziale Union)",1,0)
data$greens <- ifelse(data$party_2021_2=="GRÜNE (Bündnis 90/Die Grünen)",1,0)
data$afd <- ifelse(data$party_2021_2=="AfD (Alternative für Deutschland)",1,0)



x.vars <- c("incomequintile", "male", "I(age/10)",  "eastgermany", "whitecollarworker", "I(left_right_4/2)", "climate_policyn", "climatehumans", "afd", "fdp","linke", "spd", "cdu", "greens")


y.vars <- c("supportintervention", "supportinterventionecon","supportinterventionfrail","supportinterventionworkers")
y.varnames <- c("Belief in State Intervention", "Belief in State Intervention, Prices and Growth", "Belief in State Intervention, Sick and Elderly", "Belief in State Intervention, Workers")
y.names <- c("Belief in State Intervention", "Belief in State Intervention, Just Economics","Belief in State Intervention, Just Sick and Elderly","Belief in State Intervention, Just Workers")



for(i in 1:length(y.vars)){
  regs.overallsupport <- t(sapply(x.vars, function(x) runreg_state_man_cov(y.vars[i], x))) %>% as.data.frame()
  
  colnames(regs.overallsupport) <- c("Variable", "Coef", "se")
  
  regs.overallsupport$Coef <- regs.overallsupport$Coef %>% paste %>% as.numeric()
  regs.overallsupport$se <- regs.overallsupport$se %>% paste %>% as.numeric()
  
  regs.overallsupport$lowci <- regs.overallsupport$Coef -1.96*regs.overallsupport$se
  regs.overallsupport$highci <- regs.overallsupport$Coef +1.96*regs.overallsupport$se
  regs.overallsupport$variablename <- as.factor(c("Income (in Quintiles)", "Male", "Age (in Decades)", "East Germany", "White Collar", "Right-Wing (5-point Scale)","Support Climate Policy (5-point Scale)", "Belief in Man-Made Climate Change", "AfD Voter", "FDP Voter", "Linke Voter", "SPD Voter", "CDU Voter", "Greens Voter"))
  
  regs.overallsupport$variablename <- factor(regs.overallsupport$variablename,
                                             levels = levels(regs.overallsupport$variablename)[
                                               c(2,8,10, #demographics
                                                 5,14, #other circumstances
                                                 3,13,11, #beliefs
                                                 1,4,6,7,9,12)
                                             ])
  
  regs.list <- lapply(x.vars,function(x) runreg_state_man_table(y.vars[i], x))
  

  
  table <- modelsummary(regs.list, gof_omit = "AIC|BIC|Log.Lik|RMSE",
               coef_rename = c("I(age/10)" = "Age in Decades",
                               "I(left_right_4/2)" = "Left-Right 5 pt Scale",
                               "climate_policyn" = "Support Climate Policy",
                               "climatehumans" = "Belief in MM-CC" ),
               title=paste("Binary Regressions, DV:", y.varnames[i]),
               output="latex"
               )
  
  header.title<- paste("DV:", y.varnames[i])
  
  
  tableheader <- c(" " = 1,  header.title = 14)
  
  # set vector names 
  names(tableheader) <- c(" ", header.title)
  
  
  regtable <- table %>% kableExtra::add_header_above(header=tableheader) 


  save_kable(regtable, paste0("Tables/bivariates_coefs_state_management_",y.vars[i],".tex"))
  
  plot.coef <- ggplot(data=regs.overallsupport)+
    geom_point(aes(x=Coef, y=variablename))+
    geom_segment(aes(y=variablename , x=lowci, xend=highci,
                     yend=variablename))+
    geom_vline(aes(xintercept=0), linetype=2)+
    theme_bw()+
    xlab("Coefficient Bivariate Regression")+
    ylab("Independent Variable")+
    ggtitle(paste0("Bivariate Regression Coefficients \n DV:", y.names[i]))+
    theme(text = element_text(size=14),
          plot.title = element_text(size=18, hjust = 0.5),
          axis.text.y = element_text(size=16),
          legend.position = "bottom",
          legend.title = element_text(size=16), #change legend title font size
          legend.text = element_text(size=16),
          axis.title= element_text(size=16))
  
  print(plot.coef)
  
  ggsave(filename = paste0("Graphs/bivariates_coefs_state_management_",y.vars[i],".png"),
         plot = plot.coef,
         width = 12, height = 8)
}


####histograms Support Intervention

hist_distribution_function <- function(x, legendtext, title, binwidth){
  
d <-   data[,x] %>% as.vector
d <- d[[1]]

ggplot(data=data)+
    geom_histogram(aes(x=d), binwidth = binwidth, color="black", fill="lightgrey")+
    geom_vline(aes(xintercept = median(d,na.rm=T)), linetype="dashed", size=1.5)+
    theme_bw()+
    xlab(legendtext)+
    ylab("Count")+
    ggtitle(title)+
    theme(text = element_text(size=14),
          plot.title = element_text(size=18, hjust = 0.5),
          axis.text = element_text(size=14),
          legend.position = "bottom",
          legend.title = element_text(size=14), #change legend title font size
          legend.text = element_text(size=14),
          axis.title= element_text(size=14))
}

main <- hist_distribution_function("supportintervention", "Mean Belief in State Intervention", "Distribution of Mean Response Belief in State Intervention Questions", 0.05)
main
st1 <- hist_distribution_function("sti1", "Provide jobs for those willing to work", "", 0.5)
st1
st2 <- hist_distribution_function("sti2", "Provide an adequate standard of living for the unemployed", "", 0.5)
st2
st3 <- hist_distribution_function("sti3", "Keep prices under control", "", 0.5)
st3
st4 <- hist_distribution_function("sti4", "Support industry to ensure economic growth", "", 0.5)
st4
st5 <- hist_distribution_function("sti5", "Provide healthcare for the sick", "", 0.5)
st5
st6 <- hist_distribution_function("sti6", "Secure an appropriate standard of living for the elderly", "", 0.5)
st6

library(grid)
plot.distributions <- grid.arrange(st1, st2, st3, st4, st5, st6, nrow=3, top= textGrob("Individual Answer Distribution: The State should...", gp=gpar(fontsize=22)))

ggsave(filename = paste0("Graphs/hist_mean_belief_se.png"),
       plot = main,
       width = 12, height = 8)

ggsave(filename = paste0("Graphs/hist_individual_beliefs_se.png"),
       plot = plot.distributions,
       width = 14, height = 16)


#table on standard deviation
sd_data <- data %>% group_by(party_2021_2) %>% dplyr::summarize(`SD Belief State Intervention` = sd(supportintervention, na.rm=T))
sd(data$supportintervention, na.rm=T)

#graph of density plots for partisanship

graph_data <- data %>% filter(!is.na(party_2021_2))

graph_data$party <- recode(graph_data$party_2021_2,
                           "FDP (Freie Demokratische Partei)" = "FDP",
                           "LINKE (Die Linke)" = "The Left",
                           "SPD (Sozialdemokratische Partei Deutschlands)" = "SPD",
                           "CDU/CSU (Christlich Demokratische Union/Christlich-Soziale Union)" = "CDU",
                           "GRÜNE (Bündnis 90/Die Grünen)" = "The Greens",
                           "AfD (Alternative für Deutschland)" = "AfD",
                           "Ich werde nicht wählen" = "Non-Voter",
                           "Andere Partei" = "Other",
                           "Ich bin nicht wahlberechtigt" = "Not Eligible")

density_plots_party <- ggplot(graph_data)+
  geom_density(aes(x=supportintervention), fill = "lightgrey")+
  geom_vline(aes(xintercept = mean(supportintervention,  na.rm=T)), linewidth=2)+
  facet_wrap(~party)+
  theme_bw()+
  xlab("Belief in State Intervention Continous Measure")+
  ggtitle("Distribution of Belief in State Intervention Continous Measure \n by Party Vote 2021")+
  theme(text = element_text(size=14),
        plot.title = element_text(size=18, hjust = 0.5),
        axis.text = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=14),
        axis.title= element_text(size=14))

ggsave(filename = paste0("Graphs/density_plots_bse_party.png"),
       plot = density_plots_party,
       width = 10, height = 10)


graph_data$leftright <- ifelse(graph_data$left==1,"left", "right")

density_plots_left_right <- ggplot(graph_data)+
  geom_density(aes(x=supportintervention), fill = "lightgrey")+
  geom_vline(aes(xintercept = mean(supportintervention,  na.rm=T)), linewidth=2)+
  facet_wrap(~leftright)+
  theme_bw()+
  xlab("Belief in State Intervention Continous Measure")+
  ggtitle("Distribution of Belief in State Intervention Continous Measure \n by Left-Right Leaning")+
  theme(text = element_text(size=14),
        plot.title = element_text(size=18, hjust = 0.5),
        axis.text = element_text(size=14),
        legend.position = "bottom",
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=14),
        axis.title= element_text(size=14))

density_plots_left_right

ggsave(filename = paste0("Graphs/density_plots_bse_leftright.png"),
       plot = density_plots_left_right,
       width = 10, height = 5)
